#=============================================================================#
#
# PROJECT:        Who Pays for Peace?
# AUTHORS:        ** anonymized for review **
# CONTACT:        ** anonymized for review **
# LAST MODIFIED:  July 5, 2022
# 
#=============================================================================#
#
# This R file contains the code used to calculate the needed sample size for T2
# 
#=============================================================================#

# Initial settings ------------------------------------------------------------

#rm(list=ls())
#getwd()

## Install and load all necessary packages ------------------------------------
# ipak function: install and load multiple R packages.
# check to see if packages are installed. Install them if they are not, then load them into the R session.

ipak <- function(pkg){  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}
packages <- c("pwr", "effectsize") 
ipak(packages)


## Load the data --------------------------------------------------------------
mydata <- readRDS("data/T1-clean-data.rds")


###### A PRIORI POWER CALCULATIONS ###### 
## 1. Models and Anova Tables to get f's ------------------------------------------
fit1 <- lm(compensation ~ comp_condition, data = mydata)
anova_table1 <- anova(fit1)
effectsize(anova_table1, type = "f")

fit2 <- lm(punishment ~ punish_condition, data = mydata)
anova_table2 <- anova(fit2)
effectsize(anova_table2, type = "f")

## 2. Power calculations ----------------------------------------------------------
pwr.anova.test(n = NULL, k = 5, f = 0.67, sig.level = 0.05, power = .95) # based on full model
pwr.anova.test(n = NULL, k = 4, f = 0.79, sig.level = 0.05, power = .95) # based on full model

# Given our remarkably large effect sizes a sample size of about 50 (=10*5) 
# would even be sufficient to obtain a power of .95

effectsize(fit1) # beta (approx. cohen's d)
effectsize(fit2) # beta (approx. cohen's d)
pwr.r.test(n = NULL, r = 0.16, sig.level = 0.05, power = .80) # based on smallest, sign. beta
pwr.r.test(n = NULL, r = 0.16, sig.level = 0.05, power = .95) # based on smallest, sign. beta